clear all; clc; close all;
%Updated: 14/08/215

v0    =load('v.txt');
h0    =load('x.txt');
gn0   =load('gn.txt');
ct0   =load('ctt.txt');
vcon0 = v0(:,2);
vaut0 = v0(:,3);

q_paid0=load('q_paid.txt');
q10    =load('q.txt');
qpaid0 = q_paid0(:,2);
q1     = q10(:,1);


bounds0 = load('bounds.txt');
boundsa = bounds0(2,:);
b       = load('bgrid.txt');
a       = load('agrid.txt');


b_next0 =load('b_next.txt');
def0    = load('default.txt');

nb=size(b,1);
na=size(vcon0,1)/nb;

dif=(boundsa(2)-boundsa(1))/(na-1);
a=boundsa(1):dif:boundsa(2);

omegac = 0.31;
mu = 1;


for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((j-1)*na + i);
        vaut(i,j)  =vaut0((j-1)*na + i);
        def(i,j)   =def0((j-1)*na + i);
        qpaid(i,j) =qpaid0((j-1)*na + i);
        q(i,j)     =q1((j-1)*na + i);
        h(i,j)     =h0((j-1)*na + i);
        gn(i,j)    =gn0((j-1)*na + i);
        cn(i,j)    = 1d0-gn(i,j);
        ct(i,j)    =ct0((j-1)*na + i);
        pn(i,j)    = (1-omegac)/omegac*(ct(i,j)/cn(i,j)).^(1+mu);
        bp(i,j)    =b_next0((j-1)*na + i);
        
        if ((vcon(i,j)<=-1d2)&&(def(i,j)==0));
       % q(i,j)     =NaN;
        h(i,j)     =NaN;
        gn(i,j)    =NaN;
        cn(i,j)    =NaN;
        ct(i,j)    =NaN;
        pn(i,j)    =NaN;
        bp(i,j)    =NaN;
        end;
    end;
end;


% Unconditional mean and std dev of y:
% E(y) = e^{mu+0.5sigma^2} = 1.0010617
% stddev(y) = sqrt [(e^{sigma^2}-1)e^{sigma^2}]= 0.04614

lowa = round(na/4)+3;   %uncond mean-1std dev
meda = round(na/2);
hgha = round(3*na/4)-3; %uncond mean -1std dev


figure;plot(b,q(lowa,:),'-r',b,q(meda,:),'-b',b,q(hgha,:),'--r');


figure;contourf(b(1:nb),a(1:na),def(:,:),20);
xlabel('b');
ylabel('a');

figure;
subplot(2,2,1);
plot(b(1:nb),bp(1:2:na/2,:));
xlabel('b'); title('bprime low cT') ;
subplot(2,2,2);
plot(b(1:nb),bp(na/2:2:na,:));
xlabel('b');title('bprime low cT') ;
subplot(2,2,3);
plot(b(1:nb),pn(1:2:na/2,:));
xlabel('b');title('pn low cT') ;
subplot(2,2,4);
plot(b(1:nb),pn(na/2:2:na,:));
xlabel('b');title('pn high cT') ;

figure;
subplot(2,2,1);
plot(b(1:nb),vcon(1:2:na/2,:));
xlabel('b');title('vcon low cT') ;
subplot(2,2,2);
plot(b(1:nb),vcon(na/2:2:na,:));
xlabel('b');title('vcon high cT');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('b');title('value functions ');
% axis([b(1) b(nb) -35 -15]);
subplot(2,2,3);
plot(b(1:nb),h(1:2:na/2,:));
xlabel('b');title('h low cT') ;
subplot(2,2,4);
plot(b(1:nb),h(na/2:2:na,:));
xlabel('b');title('h high cT') ;



figure;
subplot(2,2,1);
plot(b(1:nb),gn(1:2:na/2,:));
xlabel('b');title('gn low cT') ;
subplot(2,2,2);
plot(b(1:nb),gn(na/2:2:na,:));
xlabel('b');title('gn high cT');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('b');title('value functions ');
% axis([b(1) b(nb) -35 -15]);
subplot(2,2,3);
plot(b(1:nb),ct(1:2:na/2,:));
xlabel('b');title('ct low cT') ;
subplot(2,2,4);
plot(b(1:nb),ct(na/2:2:na,:));
xlabel('b');title('ct high cT') ;
dbstop;
dbstop;

figure;
subplot(2,1,1);
plot(a(1:na),vcon(:,1:2:nb/2),a(1:na),vaut(:,1),':');
xlabel('a');
subplot(2,1,2);
plot(a(1:na),vcon(:,nb/2:2:nb),a(1:na),vaut(:,1),':');
xlabel('a');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('b');title('value functions ');
% axis([b(1) b(nb) -35 -15]);


figure;
contourf(b(1:nb),a(1:na),def(:,:),20);
xlabel('b');
ylabel('a');



figure;
subplot(2,2,1);
plot(b(1:nb),bp(1:na/2,:));
xlabel('b');
subplot(2,2,2);
plot(b(1:nb),bp(na/2:na,:));
xlabel('b');
subplot(2,2,3);
plot(b(1:nb),q(1:na/2,:));
xlabel('b');
subplot(2,2,4);
plot(b(1:nb),q(na/2:na,:));
xlabel('b');
%plot(b(1:nb),bp(lowa,:),'-r',b(1:nb),bp(meda,:),':r',b(1:nb),bp(hgha,:),'--r','LineWidth',1.5);
xlabel('b');
title('debt policy');



% a   =load('a.txt');
% b   =load('b.txt');
% 
% 
% p0  =load('gp.txt');
% h0  =load('gh.txt');
% ibp0=load('gibp.txt');
% ct0 =load('gct.txt');
% m0  =load('gm.txt');
% gn0 =load('ggn.txt');
% tau0=load('gtau.txt');
% w0  =load('gw.txt');
% def0=load('gdef1.txt');
% q0  =load('gq.txt');
% 
% paut0  =load('gpaut.txt');
% haut0  =load('ghaut.txt');
% ctaut0 =load('gctaut.txt');
% maut0  =load('gmaut.txt');
% gnaut0 =load('ggnaut.txt');
% tauaut0=load('gtauaut.txt');
% waut0  =load('gwaut.txt');
% 
% na=rows(a);
% nk=1;
% ns=na/nk;
% nb=rows(b);
% 
% alpha = 0.63;
% kappa  =0.133;
% pm     = 1;
% gamma  = 0.43;
% omegac = 0.3;
% yloss  = 0.98;
% 
% atrans = ones(nb,1)*a';
% mstar = (pm./(a.*gamma)).^(1/(gamma-1));
% 
% for i=1:na;
%     for j=1:nb;
%         vcon(i,j)  =vcon0((i-1)*nb + j);
%         
%         p(i,j)  =p0((i-1)*nb + j);
%         h(i,j)  =h0((i-1)*nb + j);
%         ibp(i,j)=ibp0((i-1)*nb + j);
%         ct(i,j) =ct0((i-1)*nb + j);
%         m(i,j)  =m0((i-1)*nb + j);
%         gn(i,j) =gn0((i-1)*nb + j);       
%         tau(i,j)=tau0((i-1)*nb + j);       
%         w(i,j)  =w0((i-1)*nb + j);   
%         def(i,j)=def0((i-1)*nb + j);   
%         q(i,j)  =q0((i-1)*nb + j);   
%         rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
%     end;
% end;
% 
% cn =h.^alpha-gn;
% const = (ct.^omegac).*(cn.^(1d0-omegac));
% bp = b(ibp);
% 
% for i=1:na;
%     vaut(i)  =vaut0(i);
%     
%     paut(i)  =paut0(i);
%     haut(i)  =haut0(i);
%     ctaut(i) =ctaut0(i);
%     maut(i)  =maut0(i);
%     gnaut(i) =gnaut0(i);       
%     tauaut(i)=tauaut0(i);       
%     waut(i)  =waut0(i);   
% end;
% 
% 
% cnaut =haut.^alpha-gnaut;
% constaut = (ctaut.^omegac).*(cnaut.^(1d0-omegac));
% 
% for i=1:na;
%     temp1 = find(ct(i,:) > 0);
%     temp2 = find(cn(i,:) > 0);
%     temp3 = find(p(i,:) > 0);
% 
%     iminct(i) = temp1(1);
%     imincn(i) = temp2(1);
%     iminp(i) = temp3(1);
%     
%     imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
%     
%     v(i,1:imin(i)-1)  = NaN; 
%     p(i,1:imin(i)-1)  = NaN; 
%     bp(i,1:imin(i)-1) = NaN; 
%     ibp(i,1:imin(i)-1)= NaN; 
%     ct(i,1:imin(i)-1) = NaN; 
%     cn(i,1:imin(i)-1) = NaN; 
%     const(i,1:imin(i)-1) = NaN;     
%     gn(i,1:imin(i)-1) = NaN; 
%     h(i,1:imin(i)-1)  = NaN; 
%     m(i,1:imin(i)-1)  = NaN; 
%     tau(i,1:imin(i)-1)= NaN; 
% 	rec(i,1:imin(i)-1)= NaN;
% %    def(i,1:imin(i)-1)= NaN; 
% 	q(i,1:imin(i)-1)  = NaN;
% 
% %     temp1 = find(ctaut(i,:) > 0);
% %     temp2 = find(cnaut(i,:) > 0);
% %     temp3 = find(paut(i,:) > 0);
% % 
% %     iminctaut(i) = temp1(1);
% %     imincnaut(i) = temp2(1);
% %     iminpaut(i)  = temp3(1);
% % 
% %     imin(i)  = max(iminctaut(i),max(imincnaut(i),iminpaut(i)));
% % 
% %     vaut(i,1:imin(i)-1)  = NaN; 
% %     paut(i,1:imin(i)-1)  = NaN; 
% %     ctaut(i,1:imin(i)-1) = NaN; 
% %     cnaut(i,1:imin(i)-1) = NaN; 
% %     constaut(i,1:imin(i)-1) = NaN;     
% %     gnaut(i,1:imin(i)-1) = NaN; 
% %     haut(i,1:imin(i)-1)  = NaN; 
% %     maut(i,1:imin(i)-1)  = NaN; 
% %     tauaut(i,1:imin(i)-1)= NaN; 
% % 
% end;
% 
% for i=1:ns;
% for j=1:nk;        
%     vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
%     p2(i,j,:)  = p((i-1)*nk+j,:);  
%     bp2(i,j,:) = bp((i-1)*nk+j,:);  
%     ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
%     ct2(i,j,:) = ct((i-1)*nk+j,:);  
%     cn2(i,j,:) = cn((i-1)*nk+j,:);  
%     const2(i,j,:) = const((i-1)*nk+j,:);      
%     gn2(i,j,:) = gn((i-1)*nk+j,:);  
%     h2(i,j,:)  = h((i-1)*nk+j,:);  
%     m2(i,j,:)  = m((i-1)*nk+j,:);  
%     tau2(i,j,:)= tau((i-1)*nk+j,:);  
% 	rec2(i,j,:)= rec((i-1)*nk+j,:); 
%     def2(i,j,:)= def((i-1)*nk+j,:);  
% 	q2(i,j,:)  = q((i-1)*nk+j,:); 
% 
%     vaut2(i,j)  = vaut((i-1)*nk+j); 
%     paut2(i,j)  = paut((i-1)*nk+j);  
%     ctaut2(i,j) = ctaut((i-1)*nk+j);  
%     cnaut2(i,j) = cnaut((i-1)*nk+j);  
%     constaut2(i,j) = constaut((i-1)*nk+j);      
%     gnaut2(i,j) = gnaut((i-1)*nk+j);  
%     haut2(i,j)  = haut((i-1)*nk+j);  
%     maut2(i,j)  = maut((i-1)*nk+j);  
%     tauaut2(i,j)= tauaut((i-1)*nk+j);  
% end;
% end;
% 
% 
% w2    = alpha*p2.*h2.^(alpha-1);
% waut2 = alpha*paut2.*haut2.^(alpha-1);
% % 
% % sigma^2 = 0.029^2/(1-0.777^2)
% % mu = 0
% 

% 
% 
% if (nk==2);
% %default decision and prices
% def21(:,:) = def2(:,1,:);
% def22(:,:) = def2(:,2,:);
% 
% price21(:,:) = q2(:,1,:);
% price22(:,:) = q2(:,2,:);
% 
% figure;
% subplot(2,2,1);
% contourf(b(1:nb),a(1:ns),def21(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{L}');
% subplot(2,2,2);
% contourf(b(1:nb),a(1:ns),def22(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{H}');
% subplot(2,2,[3 4]);
% plot(b(1:nb),price21(lowa,:),'-r',b(1:nb),price21(meda,:),':r',b(1:nb),price21(hgha,:),'--r',...
%     b(1:nb),price22(lowa,:),'-b',b(1:nb),price22(meda,:),':b',b(1:nb),price22(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');

% %%
% %value functions: repayment vs autarky
% vcon21(:,:) = vcon2(:,1,:);
% vcon22(:,:) = vcon2(:,2,:);
% vaut21(:,:) = vaut2(:,1);
% vaut22(:,:) = vaut2(:,2);
% 
% figure;
% subplot(2,1,1);
% plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
%     vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions \kappa_{L}');
% axis([b(1) b(nb) -35 -15]);
% subplot(2,1,2);
% plot(b(1:nb),vcon22(lowa,:),'-r',b(1:nb),vcon22(meda,:),':r',b(1:nb),...
%     vaut22(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut22(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions \kappa_{H}');
% axis([b(1) b(nb) -35 -15]);
% end;
% %%
% %optimal policies and prices: repayment vs autarky
% NameArray = {'LineStyle','Color','LineWidth'};
% ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';
% 
% NameArray0 = {'LineStyle','Color','LineWidth'};
% ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
% 
% [ilowa_aut]=find(def2(lowa,1,:)==0);
% [ilowa_rep]=find(def2(lowa,1,:)==1);
% h2(lowa,1,ilowa_aut) = NaN;
% ha2(lowa,1,1:nb) = haut2(lowa,1);
% ha2(lowa,1,ilowa_rep) = NaN;
% 
% [ihgha_aut]=find(def2(hgha,1,:)==0);
% [ihgha_rep]=find(def2(hgha,1,:)==1);
% h2(hgha,1,ihgha_aut) = NaN;
% ha2(hgha,1,1:nb) = haut2(hgha,1);
% ha2(hgha,1,ihgha_rep) = NaN;
% 
% p2(lowa,1,ilowa_aut) = NaN;
% pa2(lowa,1,1:nb) = paut2(lowa,1);
% pa2(lowa,1,ilowa_rep) = NaN;
% 
% p2(hgha,1,ihgha_aut) = NaN;
% pa2(hgha,1,1:nb) = paut2(hgha,1);
% pa2(hgha,1,ihgha_rep) = NaN;
% 
% gn2(lowa,1,ilowa_aut) = NaN;
% gna2(lowa,1,1:nb) = gnaut2(lowa,1);
% gna2(lowa,1,ilowa_rep) = NaN;
% 
% gn2(hgha,1,ihgha_aut) = NaN;
% gna2(hgha,1,1:nb) = gnaut2(hgha,1);
% gna2(hgha,1,ihgha_rep) = NaN;
% 
% bp2(lowa,1,ilowa_aut) = NaN;
% bp2(hgha,1,ihgha_aut) = NaN;
% 
% cn2(lowa,1,ilowa_aut) = NaN;
% cna2(lowa,1,1:nb) = cnaut2(lowa,1);
% cna2(lowa,1,ilowa_rep) = NaN;
% 
% cn2(hgha,1,ihgha_aut) = NaN;
% cna2(hgha,1,1:nb) = cnaut2(hgha,1);
% cna2(hgha,1,ihgha_rep) = NaN;
% 
% ct2(lowa,1,ilowa_aut) = NaN;
% cta2(lowa,1,1:nb) = ctaut2(lowa,1);
% cta2(lowa,1,ilowa_rep) = NaN;
% 
% ct2(hgha,1,ihgha_aut) = NaN;
% cta2(hgha,1,1:nb) = ctaut2(hgha,1);
% cta2(hgha,1,ihgha_rep) = NaN;
% 
% m2(lowa,1,ilowa_aut) = NaN;
% ma2(lowa,1,1:nb) = maut2(lowa,1);
% ma2(lowa,1,ilowa_rep) = NaN;
% 
% m2(hgha,1,ihgha_aut) = NaN;
% ma2(hgha,1,1:nb) = maut2(hgha,1);
% ma2(hgha,1,ihgha_rep) = NaN;
% 
% w2(lowa,1,ilowa_aut) = NaN;
% wa2(lowa,1,1:nb) = waut2(lowa,1);
% wa2(lowa,1,ilowa_rep) = NaN;
% 
% w2(hgha,1,ihgha_aut) = NaN;
% wa2(hgha,1,1:nb) = waut2(hgha,1);
% wa2(hgha,1,ihgha_rep) = NaN;
% 
% tau2(lowa,1,ilowa_aut) = NaN;
% taua2(lowa,1,1:nb) = tauaut2(lowa,1);
% taua2(lowa,1,ilowa_rep) = NaN;
% 
% tau2(hgha,1,ihgha_aut) = NaN;
% taua2(hgha,1,1:nb) = tauaut2(hgha,1);
% taua2(hgha,1,ihgha_rep) = NaN;
% 
% figure;
% subplot(4,2,1);P = plot(b,squeeze(h2(lowa,1,:)),b,squeeze(ha2(lowa,1,:)),b,squeeze(h2(hgha,1,:)),b,squeeze(ha2(hgha,1,:)));
% title('h');xlim([b(1) b(end)]); ylim([0.7 1.01]);% -0.05 1.05]);
% set(P,NameArray,ValueArray);
% subplot(4,2,2);P = plot(b,squeeze(p2(lowa,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(p2(hgha,1,:)),b,squeeze(pa2(hgha,1,:)));
% title('p');xlim([b(1) b(end)]);%axis([b(1) b(end) 0.495 0.7]);
% set(P,NameArray,ValueArray);
% subplot(4,2,3);P = plot(b,squeeze(gn2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gn2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
% title('g_{N}'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.05 1.05]);
% set(P,NameArray,ValueArray);
% subplot(4,2,4);P = plot(b,squeeze(bp2(lowa,1,:)),b,squeeze(bp2(hgha,1,:)));
% title('b''');xlim([b(1) b(end)]);%axis([b(1) b(end) -1.2 0]);
% set(P,NameArray0,ValueArray0);
% subplot(4,2,5);P = plot(b,squeeze(cn2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cn2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
% xlim([b(1) b(end)]);%axis([b(1) b(end) 0.45 1.05]);
% set(P,NameArray,ValueArray);
% title('c_{N}');
% subplot(4,2,6);P = plot(b,squeeze(ct2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ct2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
% title('c_{T}');xlim([b(1) b(end)]);%axis([b(1) b(end) 0.1 0.3]);
% set(P,NameArray,ValueArray);
% subplot(4,2,7);P = plot(b,squeeze(w2(lowa,1,:)),b,squeeze(wa2(lowa,1,:)),b,squeeze(w2(hgha,1,:)),b,squeeze(wa2(hgha,1,:)));
% title('w');xlim([b(1) b(end)]);%axis([b(1) b(end) -0. 0.3]);
% set(P,NameArray,ValueArray);
% subplot(4,2,8);P = plot(b,squeeze(tau2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(tau2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
% title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
% set(P,NameArray,ValueArray);
% 
% %%
% % plot time series
% seriesp   = load('seriesp.txt');
% seriesr   = load('seriesr.txt');
% seriesh   = load('seriesh.txt');
% seriesb   = load('seriesb.txt');
% seriesct  = load('seriesct.txt');
% seriesm   = load('seriesm.txt');
% seriesgn  = load('seriesgn.txt');
% seriestau = load('seriestau.txt');    
% seriesa   = load('seriesa.txt');    
% seriesbor = load('seriesbor.txt');    
% 
% n = length(seriesp);
% 
% seriescn= seriesh.^alpha-seriesgn;
% seriesy = seriesa.*seriesm.^gamma+seriesp.*seriesh.^alpha; % -pm*seriesm;
% seriesgny = seriesgn./seriesy;
% 
% time =1:n;
% 
% ValueArray1 = {'-';'blue';1.5}';
% ValueArray2 = {'-','--';'red','blue';1.5,1.5}';
% ValueArray3 = {'-','--','-';'red','blue','black';1.5,1.5,1.5}';
% 
% % No shaded areas for default episodes
% % figure;
% % subplot(3,2,1); P=plot(time,seriesa');title('a');
% % set(P,NameArray,ValueArray1);
% % subplot(3,2,2);P=plot(time,seriesp);title('p');
% % set(P,NameArray,ValueArray1);
% % subplot(3,2,3);P=plot(time,seriesct);
% % %AX=legend('c_{T}','Location','NorthWest','Orientation','Horizontal');
% % title('c_T');
% % set(P,NameArray,ValueArray1);
% % %LEG = findobj(AX,'type','text');
% % %set(LEG,'FontSize',7);
% % subplot(3,2,4); PPP=plot(time,seriescn,time,seriesh,time,seriesgn);
% % axis([1 n -0.05 1.05]);
% % AX=legend('h','c_{N}','g_{N}','Location','NorthWest','Orientation','Horizontal');title('h,c_{N} & g_{N}');
% % set(PPP,NameArray,ValueArray3);
% % LEG = findobj(AX,'type','text');
% % set(LEG,'FontSize',7);
% % subplot(3,2,5); P=plot(time,seriesr);title('spreads');
% % set(P,NameArray,ValueArray1);axis([1 n -0.01 0.1]);
% % subplot(3,2,6); P=plot(time,-seriesb);title('b');
% % set(P,NameArray,ValueArray1);
% 
% 
% seriesbor = 1-seriesbor;
% seriesbor(199)=1;
% seriesbor(200)=0;
% seriesh(200) = 0;
% ntime = 200;
% 
% figure;
% subplot(3,2,1); 
% shadedTimeSeries(time,seriesa',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('a');
% AX1 = gca;box on;
% set(AX1,'YTick',0.85:0.1:1.15);
% axis(AX1, [time(1) time(ntime) 0.85 1.15]);
% subplot(3,2,2);
% shadedTimeSeries(time,seriesp',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('p');
% AX1 = gca;box on;
% set(AX1,'YTick',0:0.2:1);
% axis(AX1, [time(1) time(ntime) 0.95*min(seriesp) 1.05*max(seriesp)]);
% subplot(3,2,3);
% shadedTimeSeries(time,seriesct',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('c_T');
% AX1 = gca;box on;
% set(AX1,'YTick',0.05:0.05:1.15);
% axis(AX1, [time(1) time(ntime) 0.9*min(seriesct)  1.1*max(seriesct)]); 
% subplot(3,2,4); 
% shadedTimeSeries(time,seriesh',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('h,c_N & g_N');
% AX1 = gca;box on;
% axis(AX1, [time(1) time(ntime) 0 1.05]);
% hold on;
% H0 = line(time,seriesh','Color','b','LineWidth',2);
% H1 = line(time,seriescn','Color','r','LineWidth',2,'LineStyle','--');
% H2 = line(time,seriesgn','Color','g','LineWidth',2,'LineStyle',':');
% h = [H0, H1,H2];
% legend(h,'h','c_N','g_N','Location','SouthEast','Orientation','Horizontal');
% AX1 = gca;box on;
% axis(AX1, [time(1) time(ntime) 0 1.05]);
% subplot(3,2,5); 
% shadedTimeSeries(time,seriesr',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('spreads');
% AX1 = gca;box on;
% set(AX1,'YTick',0.0:0.025:0.1);
% axis(AX1, [time(1) time(ntime) 0 0.1]);
% subplot(3,2,6); P=plot(time,-seriesb);title('b');
% shadedTimeSeries(time,-seriesb',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('b');
% AX1 = gca;box on;
% set(AX1,'YTick',0.0:0.02:ceil(max(-1.1*seriesb/0.001))*0.001);
% axis(AX1, [time(1) time(ntime) 0 ceil(max(-1.1*seriesb/0.001))*0.001]);